function GDOP = SingleGdop(T,S,H)
% in   目标位置T
%      测向站位置S
%      测向站高度H
% out  GDOP

sigma_theta = 1*pi/180;% 方位角测量误差
sigma_fi = 1*pi/180;% 高低角测量误差
sigma_s = 0.01; % 站址误差

x = T(1); y = T(2);
x1 = S(1);y1 = S(2);
x2 = S(3);y2 = S(4);
x3 = S(5);y3 = S(6);
% 观测站高度
z1 = H;
z2 = z1;
z3 = z1;
% 目标高度 0 km
z = 0;

% 观测误差的协方差矩阵
dV_dVT = diag([sigma_theta^2 sigma_fi^2 ...
               sigma_theta^2 sigma_fi^2 ...
               sigma_theta^2 sigma_fi^2]);
% 计算GDOP
% 站1 
theta1 = atan((x-x1)/(y-y1));
r1 = sqrt((x-x1)^2+(y-y1)^2);
fi1 = atan((z-z1)/r1);
% 站2
theta2 = atan((x-x2)/(y-y2));
r2 = sqrt((x-x2)^2+(y-y2)^2);
fi2 = atan((z-z2)/r2);
% 站3
theta3 = atan((x-x3)/(y-y3));
r3 = sqrt((x-x3)^2+(y-y3)^2);
fi3 = atan((z-z3)/r3);
% 系数矩阵
F = [(cos(theta1))^2/(y-y1) sin(theta1)*cos(theta1)/(y1-y) 0;
    (x1-x)*cos(fi1)*sin(fi1)/(r1^2) (y1-y)*cos(fi1)*sin(fi1)/(r1^2) (cos(fi1))^2/r1
    (cos(theta2))^2/(y-y2) sin(theta2)*cos(theta2)/(y2-y) 0;
    (x2-x)*cos(fi2)*sin(fi2)/(r2^2) (y2-y)*cos(fi2)*sin(fi2)/(r2^2) (cos(fi2))^2/r2
    (cos(theta3))^2/(y-y3) sin(theta3)*cos(theta3)/(y3-y) 0;
    (x3-x)*cos(fi3)*sin(fi3)/(r3^2) (y3-y)*cos(fi3)*sin(fi3)/(r3^2) (cos(fi3))^2/r3];
h1 = [(cos(theta1))^2/(y1-y) sin(theta1)*cos(theta1)/(y-y1) 0;
    (x-x1)*cos(fi1)*sin(fi1)/(r1^2) (y-y1)*cos(fi1)*sin(fi1)/(r1^2) -(cos(fi1))^2/r1];
h2 = [(cos(theta2))^2/(y2-y) sin(theta2)*cos(theta2)/(y-y2) 0;
    (x-x2)*cos(fi2)*sin(fi2)/(r2^2) (y-y2)*cos(fi2)*sin(fi2)/(r2^2) -(cos(fi2))^2/r2];
h3 = [(cos(theta3))^2/(y3-y) sin(theta3)*cos(theta3)/(y-y3) 0;
    (x-x3)*cos(fi3)*sin(fi3)/(r3^2) (y-y3)*cos(fi3)*sin(fi3)/(r3^2) -(cos(fi3))^2/r3];
% 布站误差的协方差矩阵
dXs_dXsT = sigma_s^2.*blkdiag(h1,h2,h3)*blkdiag(h1,h2,h3)';
C = (F'*F)\F';
% 定位误差的协方差矩阵
P = C*(dV_dVT+dXs_dXsT)*C';
GDOP = sqrt(trace(P));
end